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Abstract 

In this paper we have considered a flat FLRW universe with bulk viscous Zel’dovich as the cosmic 
component. Being considered the bulk viscosity as per the Eckart formalism, we have analyzed the 
evolution of the Hubble parameter and constrained the model with the Type la Supernovae data 
thus extracting the constant bulk viscous parameter and present Hubble parameter. Further we have 
analyzed the scale factor, equation of state and deceleration parameter. The model predicts the late 
time acceleration and is also compatible with the age of the universe as given by the oldest globular 
clusters. We have also studied the phase-space behavior of the model and found that a universe 
dominated by bulk viscous Zel’dovich fluid is stable. But on the inclusion of radiation component 
in addition to the Zel’dovich fluid, makes the model unstable. Hence, even though the bulk viscous 
Zel’dovich fluid dominated universe is a feasible one, the model as such failed to predict a prior 
radiation dominated phase. 


1 Introduction 

Observational data on Type-la supernovae [UlSlEliaiHlIllEl and the CMB [8l|9] has confirmed with 
sufficient accuracy that nearly seventy percent of energy of the universe is an exotic form called dark 
energy which is responsible for the current acceleration of the universe. The remaining part of the 
cosmic components consists of nearly twenty three to twenty four percentage of weakly interacting 
dark matter dni HD [la Halil] and a few percentage of luminous matter and radiation. Even with the 
overwhelming evidence for the existence of these cosmic components the current observational data 
does not rule out the possible existence of other exotic fluid components. One of the example for such 
fluids is the dark radiation, which can exist in the early or later stage during the evolution of the 
universe m- Recently, considerable attention has been paid to the study of another exotic fluid, the 
Zel’dovich fluid or stiff fluid, first studied by Zel’dovich [T6| . Zel’dovich fluid is a perfect fluid in which 
the speed of sound is equal to the speed of light so that the equation of state becomes ujz = Pz IPz = 1, 
the highest value a fluid can have in consistency with causality. 

Zel’dovich fluid or stiff fluid behavior particularly in the cosmological context had been considered 
by many. In dealing with the self interaction between dark matter components, authors in reference m 
have shown that the self interaction field will behave like a stiff fluid. The existence of Zel’dovich fluid 
was confirmed in the Horava-Lifshitz gravity based cosmological models, when the so called detailed 
balancing conditions [laHS] were relaxed [201 [2D. The relevance of the existence of the Zel’dovich 


fluid in the early universe was discussed in reference [22]. In certain inhomogeneous cosmological 
models stiff fluid was arised as an exact non singular solution |23[ 124) . In the standard evolution 
of the Friedmann universe the density of the Zel’dovich fluid is found to be decreasing faster than 
radiation and matter. So its effect on early universe would be larger. One of the phenomena that 
took place in the early universe is the primordial nucleosynthesis which might be influenced by the 
presence of stiff fluid. In reference [25] the authors have found a limit on the density of the stiff fluid 
from the constraints on the abundances of the light elements. Besides, there are no empirical facts in 
rebuttal to the stiff fluid. 

In an expanding universe there arise deviations from local thermodynamic equilibrium. Conse¬ 
quently there can arise bulk viscosity in the cosmic fluid which will restore the equilibrium |26| . This 
bulk viscosity modifies the effective pressure of the fluid in order to facilitate regaining of the equilib¬ 
rium situation. As soon as the equilibrium is reached the bulk viscous pressure vanishes |23EZj. In 
the context of inflation in the early universe it was already shown that an imperfect fluid with bulk 
viscosity can cause the early accelerated expansion |2H| . The late time viscous universe was studied in 
reference |29) . A considerable number of studies by including bulk viscosity in the dark matter setter 
were carried out by many in the context of the late acceleration of the universe laolEIlEl]. 

Recently considerable interest have been shown in the study of Zel’dovich fluid in an expanding 
universe. In reference |33j the authors studied the evolution of viscous Zel’dovich fluid in a flat universe 
and found that it can have considerable effect even in the late universe. However, they have n’t tried to 
constrain the model with cosmological model observational data to arrive at a realistic picture. In the 
present work we are trying to compare the evolution of a flat universe with bulk viscous Zel’dovich fluid 
with the latest cosmological data on Type I a supernovae. We have evaluated the model parameters 
including the transport coefficient of bulk viscosity and studied the evolution particularly in the late 
stage of the universe. The paper is organized as follows. In section 2, we derive the Hubble parameter 
and constrained it using type la Supernovae data to extract the constant bulk viscous parameter and 
the present value of the Hubble parameter. We also include the evolution of the equation of state 
parameter and deceleration in this section. In section 3, we present our analysis on the space-space 
structure of the model, followed by conclusions in section 4. 


2 The bulk viscous Zel’dovich fluid model 

The main feature of the Zel’dovich fluid is that sound velocity in the fluid is equal to that of light. 
The equation of state m is given by, 

Pz = Pz- (1) 

A similar equation of state was studied with reference to some special case by Masso and others |34j . 
In including the viscosity in the fluid we will follow the Eckart formulation which deals with the viscous 
dissipative processes occurred in a thermodynamical system when it deviates from local equilibrium. 
An equivalent formulation was developed by Landau and Lifshitz |36) . However it was noted that 
the equilibria in Eckart’s frame are unstable m and signals were propagated through the fluid at 
superluminal velocities |38j . These draw backs were rectified in a more general formalism by Israel et 
al. |39l [40] from which Eckart’s theory follows as a hrst order limit. But many authors are still using 
Eckart’s theory because of its simplicity. For example, Eckart’s formalism was used in some models 
on the late acceleration of the universe caused by the bulk viscous dark matter mmiiiis]. In the 
mean time Hiscock et al. jUj have shown that Eckart formalism can be favored over the Israel-Stewart 
formalism in inflation during the early universe using bulk viscosity. In the present work we too follow 
Eckart’s approach so that the effective pressure of the bulk viscous Zel’dovich fluid can be expressed 


as 


( 2 ) 


Pz=Pz- KH 

where Q is the coefficient of viscosity and H is the Hubble parameter. 

We consider a flat Freedmann universe with FLRW metric given as 

ds^ = dt^ — a^{t){dr^ + r^d6^ + sin^ BdcfP') (3) 

where t is the cosmic time, a{t) is the scale factor of expansion, r, 0, (j) are the comoving coordinates. 

This when combined with the Einstein’s field equations gives the dynamical equations 

= p/3 (4) 

and the conservation equation 

p + 3H{p + p) = 0 (5) 

Here we follow the standard convention SttG = 1. When these equations are combined the equation 
Q for effective pressure gives, 

dH/dt - 3H^ - = 0 (6) 

Solving these equations after changing the variable from t to x = logo we get 

ft = i {CX6 - C>-=). (7) 

where h = H/Hq, Hq is the present Hubble parameter and ^ = C/Hq is the dimensionless viscous 

parameter. For C = 0 Hubble parameter becomes H ~ HQa~^. The density of the Zel’dovich fluid 

will then evolve as p ~ and the scale factor will evolve as a ~ and hence the universe 

would have eternally decelerated and hence the effect of the Zel’dovich fluid will be relevant to the 
early epoch of the universe [25]. On the other hand, the presence of viscosity will directly contribute 
a negative term, the effect of which will depend upon the value assumed by the viscous parameter 
If C > 6 the scale factor always grows exponentially with time, in other words, an eternal acceleration. 
While for 0 < C < 6 the scale factor shows an initial deceleration followed by an acceleration in 
expansion in the later phase. So the admissible values of C is very important in this model which is 
to be evaluated by the observational constraints. 


2.1 Extraction of the model parameters using Type la supernovae 
data 


The model parameter ( and the Hubble parameter Hq can be extracted using Type la supernovae data. 
We have used Union data which consists of 307 data points |15| in the redshift range 0.01 < z < 1.6. 
The distance modulus of the supernova at a red shift z is 


p,i{z) = {m - M) = 5 logio dUz) + 25 


( 8 ) 


where m is the apparent magnitude, M is the absolute magnitude and dL{z) is the luminosity distance 
of a supernova in a flat FLRW universe, which is calculated with the relation. 


dUz) = 


c(l + z) dz 


Ha 


f 


h{z') 


(9) 
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Figure 1: Confidence intervals for the parameters ( and Hq. The outer curve corresponds to 99.99% 
probability and the inner one corresponds to 99.73 % probability. The lower dot represents the values of 
the parameters corresponding to the minimum of the 


where h{z') is identical with the normalized Hubble parameter given in equation^. The distance 
moduli of supernovae at various red shifts are calculated and are compared with the corresponding 
observational data. We then construct the statistical function 


X 


2 


n 


E 


[Pi Pfc] 


(Jt 


( 10 ) 


where is the observed distance modulus of the supernova, and ak is the variance of the mea¬ 
surement and n is the number of data points. The parameter values are obtained by minimizing the 
function. The confidence regions in figure for the parameters C and Hq are then constructed 
for 99.73% and 99.99% respectively to find the best estimate of the parameters. The values of the 
parameters are shown in table and for a comparison we have also evaluated the parameters of ACDM 


Model 

■ 

A.min 

xl^in|d.0.f. 

c 

Ho 

Bulk viscous model 

300.264 

1.011 

5.25 

70.20 

ACDM model 

300.93 

1.013 

- 

70.03 


Table 1: The best estimates of the parameters ( and Hq evaluated with the supernova type-la union data 
307 data points. But we avoided some low red shifts data, so that the net number of data used is 297. 


model. With the statistical correction the values of parameters finally become, C = 5.25 ± 0.14 and 
Ho = 70.20 ± 0.58. 
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Figure 2: Evolution of scale factor with time. Thick line corresponds to ^ = 0.002, dashed line corresponds 
to C = 5.25 (best estimate) and dotted line corresponds to ( = 6.5. 


2.2 Evolution of cosmic parameters 


The behavior of scale factor in the Zel’dovich fluid dominated universe can be obtained from the 
equation Q as 


(C-6) + 6e^(*-*°) 

C 


a{t) = 


1/3 


( 11 ) 


At sufficiently early time the scale factor can be approximated as a{t) ~ (1 + 3Ho{t — to))^) which 
implies decelerated phase, while at large times the scale factor behaves as, a{t) ~ exp(^hfo(t ~ to)) 
showing that the universe follows an accelerated epoch at a later time. The fig([^ shows evolution 
of scale factors at various choices of C- From the figure it is seen that the behavior of scale factor is 
different for ^ > 6. If one finds the expression for the age of the universe it is of the form 


to — ts 


2H, 


-1 


log 



( 12 ) 


For ^ > 6 the age is not defined and consequently the universe does not have a big bang. But for the 
cases C < 6 the universe does have a big bang. For the best estimates of the parameter the age of the 
universe is found to be around 10-12 Gy. 

The equation of state of the bulk viscous fluid can be obtained using the standard relation 




where x = Ina. Substituting for h in equation([T^ using equation^, we have 




-1 + 


2(6 - 0 
Ca3 + (6 - 0 ■ 


(13) 


(14) 


The figQ shows the variation of equation of state parameter against the redshift at various choices 
of C,. In the extreme future {z —>■ —1), u)z ~1 and hence corresponds to de Sitter universe. Otherwise 
u)z shows a strong dependence on bulk viscous parameter. For small viscosity, equation of state 
parameter remains 1 but reduces to —1 in the distant future. For ^ < 6, like Q = 5.25, the best 
estimated value, ojz is positive for z > 1. But in the subsequent evolution, it reduces to negative 
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Figure 3: Evolution of with respect to scale factor. Thick line corresponds to C = 0.002, dashed line 
corresponds to C = 5.25 and dotted line corresponds to C = 6.5 


values and finally stabilizes at —1 as z —>■ —1. This means that for z < 1 the bulk viscous Zel’dovich 
fluid mimics the quintessence nature for (^ > 6, for instance when ( = 6.5 as in the figure, when 
uJz < —1 always, which corresponds to a phantom nature. 

We have also evaluated the deceleration parameter. The basic equation of the deceleration param¬ 
eter is 


g = -l- 


H 


(15) 


Substitnting for H = H^h from equationQ we have 


q = -l + 


3(6 - 0 

Ca3 + (6 - 0 


(16) 


As seen in the figure(|^, the deceleration parameter q remains —1 for all possible Q in the distant 
future. For a small value of C the deceleration remains at 2 until a distant future when it drops down 
to —1. For the best estimated value of ^ = 5.25, the switchover from deceleration to acceleration takes 
place at about z = 0.52 which is closely in agreement with the observational constraints. The same 
will proceed with ever increasing acceleration to an asymptotic value —1 at a distant future. For 
values of C > 6.0, such as the one indicated for ^ = 6.5 as in the figure Q, q is always negative and is 
increasing as the universe expands, saturating to -1 at z —)• — 1 in the future. So the final state of the 
universe in this model is a de Sitter universe for any positive value C- 


3 Phase space perspective 

A convenient method to understand the global picture of the model is to investigate into the equivalent 
phase space. For this, first, one has to identify the phase space variables and be able to write down 
the cosmological equations as a system of autonomous differential equations. The critical points of 
these autonomous differential equations can then be correlated to the cosmological solutions. If the 
critical points were a global attractor, then the trajectories of the autonomous system constructed 
near a critical point will always be attracted towards it independent of the initial conditions. 
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Figure 4: Evolution of deceleration parameter with respect to red shift. Thick line corresponds to C = 0.002, 
dashed line corresponds to C = 5-25 and dotted line corresponds to C = 6.5 


3.1 Analysis of Zel’dovich fluid in two dimensional phase space sit¬ 
uation 

In this section the behavior of the system in the two dimensional phase space with h and as the 
coordinates is examined. The coupled differential eqnations are 

h = Ho(^-3hi}A h = (17) 


n, = -Ho [(6 {n,-i)h-c)n,-Q= Q{h, o,) (is) 

where h = H/Ho and = pz/3H^. By setting h = 0 and (iz = 0, we obtain the following three 
critical points 

(he, ^zc) = (0,1), h!,), (0.87667,1). (19) 

The first root corresponds to a static universe, while the second root depends on the instantaneous 
value of VLz and hence it is not a fixed point. The third root {hc,^zc) = (0.87667,1) corresponds to 
an expanding universe dominated by Zel’dovich fluid. If the system is stable in the neighborhood of 
a critical point, the linear perturbation in its neighborhood in phase space decays with time. The 
perturbations around the critical points must satisfy, 


e 

V. 


\dh) 0 
dQ 
dh 



( 20 ) 


Here e and r/ are perturbations in h and Vtz respectively in the neighborhood about a given critical 
point. The corresponding Jacobian is 


\ dh ) 0 
9h)o 



= Ho 


- eh^z^ 

612 ,( 1 -H,) 


-3/i2 

I +3h(l - 212,) 


( 21 ) 


where the suffix '0' implies the value at the critical point. The secular equation leads to the eigenvalues 
describing the behavior of the phase space trajectories near the equilibrium points. 

The eigenvalues corresponding to the first critical point are found to be —368.2 and 184.2. As they 
are of opposite signs the critical point is a saddle point and hence unstable. Depending on the initial 
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Figure 5: Phase space structure around the critical points. The hrst critical point (0,1) is a saddle point. 
The vector diagram clearly indicate that the trajectories approaching this point are repelled away and are 
hnally converging on to the attractor critical point on the right. 


conditions the nearby trajectories around this point may approach the saddle point, but repelled by 
it and finally approaching a possible stable attractor in the future. As in the figure the trajectories 
are turning away from the equilibrium point as and when they approach it and finally converges on 
the critical point shown on the right side of the plot. 

The second critical point is not an isolated point, but varies with As per the relationship 
between h and it represents a rectangular hyperbola with the axes h = 0 and = 0 as asymptotes. 


The eigenvalues are found to be (—184.1 


- 161.394 > 
02 


Since both the eigenvalues are negative and real. 


the neighboring trajectories will converge on to the hyperbola and hence the ’critical point’ is a stable 
one. The hyperbola along which the Vtz dependent critical point moves has the coordinate axes h = 0 
and = 0 as the asymptotes, h = —flz as the directrix and (0.9363, 0.9363) as the focus. 

The third critical point is (hc^zc) = (0.87667,1). It is observed as in figure]^ this critical point 
is a global attractor. As the attempt was made to decouple the equation 0 in the linear limit in 
the vicinity of the critical point via equation (20) it results in the eigenvalues —184.2 and 0.0. The 
resulting two eigenvalues clearly indicate the model is stable for all possible initial conditions. It 
appears that the second eigenvalue 0 is suggestive of absence of any isolated critical point and rather 
a line segment as a continuous array of critical points. However, a close examination of the vector 
field plot as in figure shows that the field directions are invariably tilted, though slightly towards the 
critical point as they approach what on low resolution seems to be a straight line towards the isolated 
critical point. It is clearly evident from the continuous plot phase space structure as shown in hgure 

So, in reality the isolated critical point exists and the '0' eigenvalue leads to a line segment as the 
best ht close to the critical point. This is clear from the fact that the straight line does not arise 
from the original procedure of setting h = 0 and = 0 without the linear approximation and rather 
results in an isolated point. So as we said earlier depending on the initial conditions the trajectories 
emanating from the surroundings of the saddle critical point are repelled away from it and they hnally 
approache the stable critical point, (0.87667,1). This, by and large implyes the stability of the universe 
dominated by the bulk viscous Zel’dovich fluid. 
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Figure 6: Vector field plot of the phase space around the third critical point (0,87667,1). The vectors 
encircled shows their continuous tilt towards the critical point. 
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Figure 7: Plot of phase space trajectories around the critical point (0,87667,1). All trajectories are con¬ 
verging to the critical point and hence it is stable. 


3.2 Analysis of Zel’dovich fluid in the three dimensional phase space 
situation 

In the case of the present model, the analysis is better realistic when we include, besides the Zel’dovich 
fluid, the conventional radiation also. The first Friedmann equation becomes 

3772 = p , + p ^ (22) 

where p^ is the radiation density. The conservation equation for the radiation component by assuming 
a pressure = p-^/3, is 

p^ + 4Hp^ = 0 ( 23 ) 


The phase-space variables are h, 0^ and 71..^ among which the third parameter is Q.y = p^/SH"^. 
The dynamical equations for these parameters are represented by the coupled differential equations 


h = P{h, Qz, ^^ 7 ) = 


(3112 -h 2n^)h 


h 


( 24 ) 








(25) 




2 ( (30^ + 2Vt^)h - ^ ) - 6/i 


^2 + C 


and 


2 ( (30^ + 2 ^ 2 ^)/i --]-2h 


^2^ — 2? (/i j ^2^;^ ^2^) — 

The critical points are obtained by setting 

/i = 0, 12;2 = 0, 12^ = 0. 




(26) 


(27) 


and they are 

0 87667 

{he, n,e, M = {—^ -, ^^ 2 ,0); (0,1,0); (1, 0.87667,0) (28) 

out of which, the first mentioned is not fixed, having h inversely proportional to the instantaneous 
value of Qz- The second critical point corresponds to static universe and third one corresponds to 
an expanding universe dominated by bulk viscous Zel’dovich fluid. It is to be noted that there is no 
critical point corresponds to a radiation dominated phase. The stability of the equilibrium points in 
the case of these three critical points are obtained (this time in the 3D phase-space case), once again 
by looking at the behavior of phase-space trajectories close to them and generated due to different 
initial conditions. The coupled differential equations in the linear limit in matrix representation, in 
the neighborhood of the equilibrium points are 


' (dP) 

\ dh ) 0 

\dh)^ 

(m\ 

\ dh ) 0 


dP 

dQz 

dQ 

dVe, 

dR 

dn,. 


dP 

d^^'y J Q 

doS 

dVl^ J g 

dR \ 

<9f!^ J g_ 


(29) 


where e, f] and h are first order perturbation terms of P{h,^}z,^'y) = h, Q{h,Qz,^'y) = and 
R{h, Q-z, ^^ 7 ) = respectively, e(t), rj{t) and u{t) being the first first order linear perturbation terms 


of h, Qz and f2.y respectively. The square matrix term in the equation (29) is the Jacobian evaluated at 


the critical point. We then decouple the differential equations (29) by means of the secular equation 
and in the process the eigenvalues corresponding to the equilibrium point ^zc, 0 ) are obtained 


as 


A 

' f ^2 

A 


-368.2 365.5712^ - ^135571 - 26920612^ -k 13364122 

-368.2 365.5712^ -h ^135571 - 26920612^ -h 13364122 


and 


A 3 — 


245.466 

12 . ’ 


(30) 

(31) 

(32) 


all depending on the instantaneous value of 12^. The critical point in this case drifts with the variation 
of 12 ^ along a rectangular hyperbola on the /i — 12 ^ plane; with the details of the hyperbola same as 
in the case of the second critical point in the section 3.1 The eigenvalues indicate the phase space 


trajectories corresponding to various initial conditions move away from the critical point and hence 
no stable situation. Even the case where 12^ = 0 the eigenvalues are such that, Ai is negative, A 2 = 0 
and A 3 positive, and so stable solution solution is implied. 

The second critical point (0,1,0) has the eigenvalues —368.2, 368.2 and 184.1 and the third critical 
point (1,0.87667, 0) has the same to be —403.778, 280.0 and 167.878. There is one negative eigenvalue 
and there are two positive eigenvalues for each critical points which again means there is no stability 
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Figure 8: Vector field plot of the phase-space structure around the critical point (1,0.87667,0) 


for the equilibrium points. This means the phase space trajectories are not attracted by any of 
the critical points in the three dimensional case. For example the vector field plot as in the figure 
clearly indicates how the phase-space trajectories corresponding to various initial conditions are 
repelled away, rather than being attracted to the second critical point. So none of the critical points 
in this case corresponds to a radiation dominated phase and even the existing critical points are not 
stable also. In fact the third critical point which corresponds to a Zel’dovich fluid dominated one, 
since it is unstable, it can be concluded that the inclusion of the radiation component may lead to a 
complete break down of the model. The bulk viscous coefficient is taken as a constant in the present 
study. Since it is a transport coefficient it may depend on the velocity of the fluid component also. 
Such a velocity dependent bulk viscous coefficient may be checked for consistency of a prior radiation 
dominated phase and that we reserve for a further work. 


4 Conclusion 

In this paper we have considered a flat universe consisting of bulk viscous Zel’dovich fluid. The 
viscosity parameter was incorporated as per the Eckart’s formalism. We have evaluated the evolution 
of the Hubble parameter. The model was constrained with SNe la data to evaluate the bulk viscous 
coefficient ( = 5.25 ± 0.14 the present value of Hubble parameter Hq = 70.20 ± 0.58. The behavior 
of the resulting scale factor shows that the model predicts a late acceleration in the expansion of the 
universe. Hence the bulk viscous Zel’dovich fluid can mimic the role of the conventional dark energy. 

We also studied the model to analyze the stability of the solutions corresponding to various sce¬ 
narios using the phase space analysis method. We first analyzed the two dimensional phase-space 


behavior, where the contribution due to radiation is neglected and found that there is a past unstable 
saddle critical point corresponding to a static universe. The phase-space trajectories originating form 
the vicinity of this saddle like point are repelled away from it and are proceeded towards the stable 
critical point corresponding to an expanding universe dominated by Zel’dovich fluid. 

In the second instance we considered a three dimensional phase-space case by incorporating the 
radiation component too. In this case no critical points are found corresponding to a prior radiation 
dominated phase and more over none of the existing critical points are stable. Hence the present 
model of the universe with bulk viscous Zel’dovich, in which bulk viscosity is characterized by a 
constant coefficient, first of all failed to predict a prior radiation dominated phase and secondly the 
very inclusion of the radiation makes the very model unstable. 
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